subroutine grid_nano
  use rbf_type
  use variable, only:ddx,ddy,ddz,ry,rz,beta,nip,Nx,Ny,Nz,dddx,dddy,dddz,nodes,xmax,ymax,zmax,pi
  use variable, only: confine,inside,nee_cube,xm,nematic_l,surface,psize,nee_cube,npatch,inside,surface,npatch
  use patch_module, only : inside_nano,surface_nano,pm,type_npatch
  use ptrajectory
  implicit none
   real(dp)      :: x,y,z,r,checkx,checky,checkz,check1,check2,R_nano
   real(dp), dimension(3) :: xs
   integer (i4b)  :: i,j,k,part


      if (npatch.eq.0) psize = 0.125
      R_nano = psize
      allocate( inside_nano(nodes))
      allocate( surface_nano(nodes))  
      allocate (pm(3*npatch) )
      allocate (type_npatch(nodes))

      inside_nano =.false.
      surface_nano =.false.

      if (npatch.ne.0)then
      write(*,160)'   Diameter nano (nm)  =',2.d0*psize*nematic_l

160  format (a,f12.4)

! Particle trajectory

   open(111,file='distance.in',status='old')
   read(111,nml = input2)
   close(111)


            
   particle:do part = 1,npatch

!      theta(part) = theta(part) *pi/180.d0 
!      phi  (part) = phi (part)  *pi/180.d0


!      pm(3*part-2) = rr(part)*sin(theta(part))*cos(phi(part))*beta   
!      pm(3*part-1) = rr(part)*sin(theta(part))*sin(phi(part))*beta   
!      pm(3*part-0) = rr(part)*cos(theta(part))               *beta   

   

      pm(3*part-2) = cmx_part(part)    *beta*2.d0    - 0.5d0*xmax
      pm(3*part-1) = cmy_part(part) *ry*beta*2.d0    - 0.5d0*ymax
      pm(3*part-0) = cmz_part(part) *rz*beta*2.d0    - 0.5d0*zmax

   write(*,*) pm(3*part-2),pm(3*part-1),pm(3*part-0)

   xs(1) = pm(3*part-2) + 0.5d0*xmax
   xs(2) = pm(3*part-1) + 0.5d0*ymax
   xs(3) = pm(3*part-0) + 0.5d0*zmax

!   xs(1) = pm(3*part-2) 
!   xs(2) = pm(3*part-1) 
!   xs(3) = pm(3*part-0) 


   do i=1,Nx
      do j=1,Ny
         do k=1,Nz
              
            if(inside(i,j,k).or.surface(i,j,k))then
               x = ddx*(i-1)
               y = ddy*(j-1)
               z = ddz*(k-1)
           
               checkx = (x-xs(1))/(R_nano-ddx/2.d0);
               checkx = checkx*checkx;
               checky = (y-xs(2))/(R_nano-ddy/2.d0);
               checky = checky*checky;
               checkz = (z-xs(3))/(R_nano-ddz/2.d0);
               checkz = checkz*checkz;
               check1 = checkx+checky+checkz;

               checkx = (x-xs(1))/(R_nano+ddx/2.d0);
               checkx = checkx*checkx;
               checky = (y-xs(2))/(R_nano+ddy/2.d0);
               checky = checky*checky;
               checkz = (z-xs(3))/(R_nano+ddz/2.d0);
               checkz = checkz*checkz;
               check2 = checkx+checky+checkz;
           
              
               if  (check1 .ge. 1.d0  .and.check2.lt.1.d0 )then
                  
                  surface_nano(nee_cube(i,j,k))=.true.
                  type_npatch(nee_cube(i,j,k))=part
               endif
         endif

        enddo
     enddo
  enddo


     do i=1,Nx
        do j=1,Ny
           do k=1,Nz

              if(inside(i,j,k).or.surface(i,j,k))then

                 x = ddx*(i-1)
                 y = ddy*(j-1)
                 z = ddz*(k-1)
           
                 checkx = (x-xs(1))/(R_nano-ddx/2.d0);
                 checkx = checkx*checkx;
                 checky = (y-xs(2))/(R_nano-ddy/2.d0);
                 checky = checky*checky;
                 checkz = (z-xs(3))/(R_nano-ddz/2.d0);
                 checkz = checkz*checkz;
                 check1 = checkx+checky+checkz;

                 checkx = (x-xs(1))/(R_nano+ddx/2.d0);
                 checkx = checkx*checkx;
                 checky = (y-xs(2))/(R_nano+ddy/2.d0);
                 checky = checky*checky;
                 checkz = (z-xs(3))/(R_nano+ddz/2.d0);
                 checkz = checkz*checkz;
                 check2 = checkx+checky+checkz;
           
              
                 if  (check1 .lt. 1.d0   )then

                    inside_nano(nee_cube(i,j,k))=.true.
 
                 endif
              endif
        enddo
     enddo
  enddo
enddo particle
endif

end subroutine grid_nano
!--------------------------------
!--------------------------------
!--------------------------------

subroutine nano_normal
  use rbf_type
  use variable, only: nn, xm, normal,planar,beta,pi,nodes,ddx,ddy,ddz,xmax,ymax,zmax,npatch
  use patch_module, only :surface_nano,pm,type_npatch
  implicit none
  integer(i4b) :: i, j, ii, jj,part
  real(dp)     :: r,zz,rr,phi,theta
  real(dp), dimension(3):: dx,dx_norm,xs
    
  
  
!  normal = 0.d0
 
  do part=1,npatch

     xs(1) = pm(3*part-2) 
     xs(2) = pm(3*part-1) 
     xs(3) = pm(3*part-0)


     do i = 1, nodes

        if (surface_nano(i))then
    
           if (type_npatch(i).eq.part)then
           ii = 3*i   


           dx(1) = xm(ii-2) -xs(1)
           dx(2) = xm(ii-1) -xs(2)
           dx(3) = xm(ii  ) -xs(3)

           
           r  = sqrt( dot_product(dx,dx) ) 
           normal(ii - 2:ii) = dx/r
             endif
        endif
     end do
  enddo
 
end subroutine nano_normal
!-----------
!-----------
!-----------

subroutine q_nano
  use rbf_type
  use variable, only : Sinitial,nodes,a,xm,direction,pi,Sbulk,planar,normal,chi,typ_part,npatch,frac_part
  use patch_module, only : surface_nano,pm,type_npatch
  use  ptrajectory, only : theta,phi
  use random_number2
  use ranftest

  implicit none
!  logical(lgt) flag
  integer(i4b) i,part,ii
  real(dp), dimension(3):: dx_norm,dx_plan
  real(dp), dimension(:,:) ,allocatable :: q
 real(dp), dimension(3):: nrandom,dx_part,norm_part,paral_part,xs,dx
    real(dp) :: r,theta_cm,phi_cm,rho_cm,dot_pr,norma

  
  allocate(q(5,nodes))


  ! Boundary nodes and Prefered anchoring

  q=0.d0 
  npatchpart:do part=1,npatch

     theta(part) = theta(part) *pi/180.d0 
     phi  (part) = phi (part)  *pi/180.d0


  do i = 1, nodes
     if(surface_nano(i))then
        if (type_npatch(i).eq.part)then
        dx_norm = normal(3*i - 2:3*i)
  !  nrandom (1) =0.d0 ; nrandom (2) = 0.d0 ; nrandom (3) = 1.d0

        nrandom (1) = sin(theta(part))*cos(phi(part))
        nrandom (2) = sin(theta(part))*sin(phi(part))
        nrandom (3) = cos(theta(part)) 

        nrandom = nrandom - normal(3*i - 2:3*i)*dot_product( nrandom, normal(3*i - 2:3*i) )

        if (nrandom(1).eq.0.d0.and.nrandom(2).eq.0.d0.and. nrandom(3).eq.0.d0)then
           dx_plan (1) = 1.d0 ;  dx_plan (2) = 0.d0 ;  dx_plan (3) = 0.d0
        else
           r    = sqrt( dot_product( nrandom, nrandom ) )
           nrandom = nrandom/r
           dx_plan = nrandom
        endif


        if(typ_part.eq.1)then  ! Pure homeotropic 
          
           call director_tensor( Sbulk, dx_norm, q(:,i) )
           call transformation( q(:,i), a(:,i) )

        else if(typ_part.eq.2)then   ! Pure Bipolar 

           call director_tensor( Sbulk, dx_plan, q(:,i) )
           call transformation( q(:,i), a(:,i) )

        else if(typ_part.eq.3)then   ! Janus Particle

! Definining the vector in spherical coordinates to control the orientation of the nanoparticle, when this is Janus type
          ! write(*,*)part,theta(part),phi(part)
           
           
           norm_part(1) = sin(theta(part))*cos(phi(part))   
           norm_part(2) = sin(theta(part))*sin(phi(part))   
           norm_part(3) = cos(theta(part)) 

              
           ! write(*,*)theta,phi
           ! write(*,*)norm_part


            !paral_part(1) = norm_part(2)
            !paral_part(2) =- norm_part(1)
            !paral_part(3) = 0.d0

            !r    = sqrt( dot_product( paral_part, paral_part ) )
            !paral_part = paral_part/r

            xs(1) = pm(3*part-2) 
            xs(2) = pm(3*part-1) 
            xs(3) = pm(3*part-0)
            
            
            ii = 3*i 
            

            if(xs(1).ne.0.d0)then
               if(xs(1).lt.0.0d0)then
                  dx(1) = xm(3*i-2) -xs(1)
               else
                  dx(1) = xm(3*i-2) -xs(1)
               endif
            else
               dx(1) = xm(3*i-2) 
            endif

            if(xs(2).ne.0.d0)then
               if(xs(2).lt.0.0d0)then
                  dx(2) = xm(3*i-1) -xs(2)
               else
                  dx(2) = xm(3*i-1) -xs(2)
               endif
            else
               dx(2) = xm(3*i-1)
            endif

            if(xs(3).ne.0.d0)then
               if(xs(3).lt.0.0d0)then
                  dx(3) = xm(3*i-0) -xs(3)
               else
                  dx(3) = xm(3*i-0) -xs(3)
               endif
            else
               dx(3) = xm(3*i-0)
            endif

            

         !   write(*,*)dx(1),xm(3*i-2),pm(3*part-2)
         !   write(*,*)dx(2),xm(3*i-1),pm(3*part-1)
         !   write(*,*)dx(3),xm(3*i-0),pm(3*part-0)


            rho_cm   = sqrt(dx(1)**2 + dx(2)**2+ dx(3)**2)
            theta_cm = acos(dx(3)/rho_cm)
            phi_cm   = atan2(dx(2),dx(1))

            norma = sqrt(dot_product(dx,dx))
!            write(*,*)norma,dx

            dx = dx / norma
            
            dot_pr = dot_product( norm_part,dx )

            !write(*,*)dot_pr,dx
!            if(theta_cm.le.pi*0.5d0)then

            if(dot_pr.gt.0.d0)then 
              dx_part = dx_plan
           else
              dx_part = dx_norm
           endif


           call director_tensor( Sbulk, dx_part, q(:,i) )
           call transformation( q(:,i), a(:,i) )

           else if(typ_part.eq.6)then   ! Ying-Yang 

! Definining the vector in spherical coordinates to control the orientation of the nanoparticle, when this is Janus type

           
            theta(part) = theta(part) *pi/180.d0 
            phi  (part) = phi (part)  *pi/180.d0



            norm_part (1) =sin(theta(part))*cos(phi(part))   
            norm_part(2) = sin(theta(part))*sin(phi(part))   
            norm_part(3) = cos(theta(part))               


            paral_part(1) = norm_part(2)
            paral_part(2) =- norm_part(1)
            paral_part(3) = 0.d0

            r    = sqrt( dot_product( paral_part, paral_part ) )
            paral_part = paral_part/r

            xs(1) = pm(3*part-2) 
            xs(2) = pm(3*part-1) 
            xs(3) = pm(3*part-0)

            

!            ii = 3*i   
            
            if(xm(3*i-2).lt.0.d0)then
               dx(1) = xm(3*i-2) +xs(1)
            else
               dx(1) = xm(3*i-2) -xs(1)
            endif

            if(xm(3*i-1).lt.0.d0)then
               dx(2) = xm(3*i-1) +xs(2)
            else
               dx(2) = xm(3*i-1) -xs(2)
            endif

            if(xm(3*i-0).lt.0.d0)then
               dx(3) = xm(3*i-0) +xs(3)
            else
               dx(3) = xm(3*i-0) -xs(3)
            endif


   !         dx(2) = xm(3*i-1) +xs(2)
   !         dx(3) = xm(3*i  ) +xs(3)

   !         dx(1) = xm(3*i-2) 
   !         dx(2) = xm(3*i-1) 
   !         dx(3) = xm(3*i  ) 

     !       write(*,*)dx

            rho_cm   = sqrt(dx(1)**2 + dx(2)**2+ dx(3)**2)
            theta_cm = acos(dx(3)/rho_cm)
            phi_cm   =atan2(dx(2),dx(1))



!!!            if(theta_cm.le.pi*0.333d0)then
!            if(theta_cm.le.(pi*(0.5d0+frac_part*cos(2.d0*phi_cm))))then

           if(part.eq.1) write(47,*)part,frac_part,phi_cm
           if(part.eq.2) write(48,*)part,frac_part,phi_cm

            if(theta_cm.le.(pi*(0.5d0+frac_part*cos(2.d0*phi_cm))))then
              dx_part = dx_plan
           else
              dx_part = dx_norm
           endif

           call director_tensor( Sbulk, dx_part, q(:,i) )
           call transformation( q(:,i), a(:,i) )

        endif

     endif
  endif
  end do

enddo npatchpart

  deallocate( q )
end subroutine q_nano
